%% Code for generating compensation-performance and flow-performance plots in optimal contract

%% Params
day_incr = 1;

dt = day_incr/365; %in days
YRS = 1;

T = round(YRS/dt);

dp0 = 0.05;
dR01 = 0.01;

y0 = 0;
p0 = round([0:dp0:1]*100)/100;
R01 = round([-1.5:dR01:1.5]*100)/100;

cases = combvec(p0, R01)'; %prior in first column and targer return in second column

n_sim = length(cases);

tt = [0:T];
Time = repmat(tt,n_sim,1);


%% Pick model
model_pick = 'lc';

V = eval(['V_',model_pick]);
Eta = eval(['Eta_',model_pick]);
Beta = eval(['Beta_',model_pick]);
C = eval(['C_',model_pick]);
K = eval(['K_',model_pick]);

mu_v = eval(['mu_v_',model_pick]);
vol_phi = eval(['vol_phi_',model_pick]);
ExRet = eval(['ExRet_',model_pick]);


%% INTERPOLANTS
V_FUN = griddedInterpolant(PHI,V);

C_FUN = griddedInterpolant(PHI,C);
K_FUN = griddedInterpolant(PHI,K);
Beta_FUN = griddedInterpolant(PHI,Beta);
Eta_FUN = griddedInterpolant(PHI,Eta);

ExRet_FUN = griddedInterpolant(PHI,ExRet);

mu_v_FUN = griddedInterpolant(PHI,mu_v);

mu_y_FUN = griddedInterpolant(PHI,mu_y);
vol_y_FUN = griddedInterpolant(PHI,vol_y);

vol_phi_FUN = griddedInterpolant(PHI,vol_phi);

%% GENERATE DATA
phi = NaN(n_sim, T+1);
phi(:,1) = cases(:,1);

dW = NaN(n_sim,T+1);

log_v = NaN(n_sim, T+1);
log_v(:,1) = 0;

c = NaN(n_sim,T+1);
k = NaN(n_sim,T+1);
eta = NaN(n_sim,T+1);
b = NaN(n_sim,T+1);
ff = NaN(n_sim,T+1);
mv = NaN(n_sim,T+1);

sphi = NaN(n_sim,T+1);

exret = NaN(n_sim,T+1);

TotRet = cases(:,2); %total return above initial expectation
dR = TotRet/T;

phi_start = phi(:,1);

for t = 1:T
    
    % Controls
    c(:,t) = C_FUN(phi_start);
    k(:,t) = K_FUN(phi_start);
    eta(:,t) = Eta_FUN(phi_start);    
    b(:,t) = Beta_FUN(phi_start);
        
    % Returns
    exret(:,t) = ExRet_FUN(phi_start);
    
    % State variables
    mv(:,t) = mu_v_FUN(phi_start);            
    sphi(:,t) = vol_phi_FUN(phi_start);
    
    %Shocks 
    dW(:,t) = (dR - exret(:,t)*dt)./sigma;
    
    %update state variables  
    phi(:,t+1) = phi_start + sphi(:,t).*dW(:,t);
    
    log_v(:,t+1) = log_v(:,t) + (mv(:,t) - 0.5*b(:,t).^2)*dt + b(:,t).*dW(:,t);
    
    phi_start = phi(:,t+1);
end

DK = log((k(:,T)./k(:,1)).*exp(log_v(:,T)));
DC = log((c(:,T)./c(:,1)).*exp(log_v(:,T)));

%% PLOTS

paperposition = [0 0 8 4.5];
definition = '-r500';
axis_font = 6;
label_font = 10;
legend_font = 8;

edgecolor = 'k';

linewidth = 1;

color1 = 'b';
color2 = '--r';

% Starting prior
ppp = 0.5;

%% Capital
VAR = DK(cases(:,1) == ppp,:) - DK(cases(:,1) == ppp & cases(:,2) == 0,:);

figure;
hold on
plot(R01, VAR, color1, 'LineWidth', linewidth);
set(gca,'FontSize',axis_font);
Xla=xlabel('Cumulative return','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$\Delta \log K$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
print('-dpng',gcf,['convexity_aum_lc.png'],'-r500');
hold off; clear h; close;





%% Compensation
VAR = DC(cases(:,1) == ppp,:) - DC(cases(:,1) == ppp & cases(:,2) == 0,:);

figure;
hold on
l1 = plot(R01, VAR, color1, 'LineWidth', linewidth);
set(gca,'FontSize',axis_font);
Xla=xlabel('Cumulative return','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$\Delta \log C$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
print('-dpng',gcf,['convexity_comp_lc.png'],'-r500');
hold off; clear h; close;
